Based on "Optimal Withdrawal Strategies for Retirement Income Portfolios" (2012) by Blanchett, Kowara, and Chen
In [150]:
%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import itertools
import functools
import operator
sns.set(color_codes=True)
In [151]:
# based on http://www.ssa.gov/OACT/STATS/table4c6.html Period Life Table, 2009
# https://drive.google.com/open?id=0Bx-R9jPONgMgX0ktTUVFbVU0TEU
lifetable_for_women = [
0.010298, # 65
0.011281,
0.012370,
0.013572,
0.014908,
0.016440, # 70
0.018162,
0.020019,
0.022003,
0.024173,
0.026706, # 75
0.029603,
0.032718,
0.036034,
0.039683,
0.043899, # 80
0.048807,
0.054374,
0.060661,
0.067751,
0.075729, # 85
0.084673,
0.094645,
0.105694,
0.117853,
0.131146, # 90
0.145585,
0.161175,
0.177910,
0.195774,
0.213849, # 95
0.231865,
0.249525,
0.266514,
0.282504,
0.299455, # 100
0.317422,
0.336467,
0.356655,
0.378055,
0.400738, # 105
0.424782,
0.450269,
0.477285,
0.505922,
0.536278, # 110
0.568454,
0.602561,
0.638715,
0.677038,
0.717660, # 115
0.760720,
0.806363,
0.851378,
0.893947 # 119
]
In [152]:
import random
def die(age):
if age >= len(lifetable_for_women): return age
if random.random() <= lifetable_for_women[age]:
return age
else:
return die(age + 1)
In [153]:
def gen_age():
return die(0) + 65
In [204]:
# Just to show what the overall mortality curve looks like
ages = []
for _ in range(250000):
ages.append(gen_age())
ages = pd.Series(ages, name='Age at death')
sns.distplot(ages)
Out[204]:
In [155]:
# import US based data from 1871 from the simba backtesting spreadsheet found on bogleheads.org
# https://www.bogleheads.org/forum/viewtopic.php?p=2753812#p2753812
dataframe = pd.read_csv('returns.csv')
dataframe.head()
Out[155]:
In [156]:
def shuffle(df, n=1, axis=0):
index = list(df.index)
random.shuffle(index)
df = df.ix[index]
df.reset_index()
return df
Generate monte carlo returns based on shuffling historical data.
The actual WER calculated is extremely sensitive to the model used for simulating returns. The more generous the model, the higher the SSR will be (since it can calculate a maximal strategy) and the lower the WER of most strategies.
This means that it is hard to compare WER's across research unless are using the same monte carlo simulations, as far as I can tell.
The relative performance of strategies should be the same across research but the absolute WER will vary.
In [206]:
def personal_returns(dataframe, age=None, inflation=False, adjustment=0.0):
df = shuffle(dataframe.copy())
if age == None:
age = gen_age()
lifetime_returns = []
for _ in range(age - 64):
row = df.iloc[_]
bonds = row['Bonds']
stocks = row['Stocks']
annual_returns = (bonds * .4) + (stocks * .6)
# The Morningstar paper isn't clear whether they use real returns or not.
# I'll assume not and leave out the inflation adjustment.
if inflation:
annual_returns -= row['Inflation']
if adjustment:
annual_returns += adjustment
lifetime_returns.append(annual_returns / 100)
return lifetime_returns
In [207]:
def prod(x):
return functools.reduce(operator.mul, x, 1)
In [208]:
def ssr(r):
def ssr_seq(r):
if len(r) == 0: return 1
else:
denom = prod(map(lambda x: x + 1, r))
return (1 / denom) + ssr_seq(r[1:])
# skip the final year of returns (since you have withdrawn the entire portfolio at the beginning of the year,
# knowing that you will die this year)
r_r = list(reversed(r[:-1]))
return 1 / ssr_seq(r_r)
In [209]:
# test that I've implemented SSR correctly.
pr = [0.067040000000000002, -0.049419999999999999, 0.20482000000000003, 0.28746000000000005]
print('SSR %f' % ssr(pr))
a0 = pr[0] + 1
a1 = pr[1] + 1
a2 = pr[2] + 1
a3 = pr[3] + 1
hand_ssr = (1.0
/ (1
+ (1/ (a0))
+ (1/ (a0 * a1))
+ (1/ (a0 * a1 * a2))
)
)
print('Hand %f' % hand_ssr)
In [210]:
def cew(cashflows):
# the paper says that results aren't sensitive to changing this (which turns out to be true!)
gamma = 4.0
def sigma(c):
return pow(c, -gamma) / gamma
constant_factor = 1.0 / len(cashflows) * gamma
base = constant_factor * sum(map(sigma, cashflows))
return pow(base, -1/gamma)
In [211]:
# test CEW
print(cew([0.05, 0.055, 0.06]))
print('0.05424659419') # calculated from a google spreadsheet
In [212]:
# starting at age 65, based on 60% domestic stocks, 40% domestic bonds
vpw_rates = map(lambda x: x/100.0, [5.0,
5.1,
5.1,
5.2,
5.3,
5.4,
5.5,
5.6,
5.7,
5.9,
6.0,
6.2,
6.3,
6.5,
6.7,
6.9,
7.2,
7.5,
7.8,
8.1,
8.5,
9.0,
9.5,
10.1,
10.9,
11.7,
12.8,
14.2,
15.9,
18.2,
21.5,
# diverge from the actual VPW...never exhaust the porfolio, just take 25% forever
25.0, #25.4,
25.0, #34.6,
25.0, #50.9,
25.0, #100
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
25.0,
])
cew(vpw_rates)
Out[212]:
In [213]:
def vpw_wer(p, s):
lifespan = len(p)
v = vpw_rates[0:lifespan]
w = cew(v) / s
return w
In [214]:
def endow_wer(p, s):
''' use a constant 5.5% withdrawal '''
lifespan = len(p)
r = list(itertools.repeat(.055, lifespan))
w = cew(r) / s
return w
In [215]:
life_expectancy = [
# from the Annuity 2000 table A-1
# starting age 65
23.02,
22.16,
21.31,
20.47,
19.63,
18.81, # 70
17.99,
17.17,
16.40,
15.62,
14.86, # 75
14.12,
13.40,
12.69,
12.01,
11.34, #80
10.70,
10.08,
9.48,
8.91,
8.37, # 85
7.85,
7.37,
6.91,
6.48,
6.08, # 90
5.72,
5.38,
5.06,
4.77,
4.50, # 95
4.24,
4.0,
3.76,
3.52,
3.29, # 100
3.05,
2.81,
2.58,
2.35,
2.12, # 105
1.91,
1.70,
1.51,
1.33,
1.16, # 110
# now just pad it out
1,
1,
1,
1,
1,
1,
1,
1,
]
rmd_rates = map(lambda x: 1.0 / x, life_expectancy)
In [216]:
def rmd_wer(p, s):
lifespan = len(p)
v = rmd_rates[0:lifespan]
w = cew(v) / s
return w
In [223]:
vpw_wers = []
endow_wers = []
rmd_wers = []
for _ in range(25000):
p = personal_returns(dataframe.copy())
s = ssr(p)
vpw_wers.append(vpw_wer(p, s))
endow_wers.append(endow_wer(p, s))
rmd_wers.append(rmd_wer(p, s))
series = pd.Series({
'VPW' : np.average(vpw_wers),
'Endowment' : np.average(endow_wers),
'RMD': np.average(rmd_wers),
})
series.plot(kind='bar')
series.head()
Out[223]:
In [222]:
vpw1_wers = []
for _ in range(25000):
p = personal_returns(dataframe.copy())
s = ssr(p)
vpw1_wers.append(vpw_wer(p, s))
vpw2_wers = []
for _ in range(25000):
p = personal_returns(dataframe.copy(), inflation=True)
s = ssr(p)
vpw2_wers.append(vpw_wer(p, s))
series = pd.Series({
'Base' : np.average(vpw1_wers),
'Inflation-adjusted' : np.average(vpw2_wers),
})
series.plot(kind='bar')
series.head()
Out[222]:
In [187]:
# chance of dying in N years of retirement
years_of_retirement = 10
1 - prod(map(lambda x: 1-x, lifetable_for_women[0:years_of_retirement]))
Out[187]:
In [188]:
pr = personal_returns(dataframe, age=65 - 1 + years_of_retirement)
s = ssr(pr)
print('Years of retirement %s. SSR %s' % (len(pr), s))
PORTFOLIO_START = 1000000
ssr_port = PORTFOLIO_START
ssr_wd = ssr_port * s
ssr_vals = []
vpw_port = PORTFOLIO_START
vpw_vals = []
for _ in range(len(pr)):
# make the withdrawal at the start of the year
ssr_port -= ssr_wd
# then add in all of that year's gains
ssr_port += ssr_port * pr[_]
ssr_vals.append(ssr_port)
vpw_wd = (vpw_rates[_] * vpw_port)
vpw_port -= vpw_wd
vpw_port += vpw_port * pr[_]
vpw_vals.append(vpw_port)
plt.plot(ssr_vals)
plt.plot(vpw_vals)
Out[188]:
In [199]:
sns.tsplot(vpw_rates[0:30], color='green')
sns.tsplot(rmd_rates[0:30], color='blue')
Out[199]:
In [203]:
vpw_vs_rmd = []
ssr_vs_rmd = []
ssrate = ssr(personal_returns(dataframe.copy(), age=95))
print('Sustainable rate for 30 years: %f' % ssrate)
for _ in range(20):
vpw_w = vpw_rates[_]
rmd_w = rmd_rates[_]
delta = vpw_w - rmd_w
vpw_vs_rmd.append(delta)
delta = ssrate - rmd_w
ssr_vs_rmd.append(delta)
sns.tsplot(vpw_vs_rmd, color='blue')
sns.tsplot(ssr_vs_rmd, color='green')
Out[203]: